rm(list=ls())

source('param_general.R')
source('param_switch.R')

fileTyp <- 'all'

print(paste0(datDir,'chambDat_minPpl_',fileTyp,'_',minZipHospPpl,'.RData'))
load(paste0(datDir,'chambDat_minPpl_',fileTyp,'_',minZipHospPpl,'.RData', sep = '')) 

chambDat[,':='(time_diff_firstBirthHosp = time_current_2_1 - time_current_1_1,
               time_diff_secBirthHosp = time_current_2_0 - time_current_1_0)]
chambDat[,time_diff_secBirth :=(time_current_2_1-time_current_2_0)]
chambDat[,time_diff_firstBirth :=(time_current_1_0-time_current_1_1)]
chambDat[,one := 1]

honoreDatUse <- subset(chambDat, age06 <= ageCut06 & noMove3==1 & hosp_id_alt_1_0 != hosp_id_alt_1_1 )
honoreDatUse[,secondHospMinFirstHosp3 := (secondHosp3==1)-(firstHosp3==1)]

##Different data sets for estimation
assignDatFunc <- function(datUse){
    datList <- list(Main = as.data.frame(datUse),
                    Normal_LD= as.data.frame(datUse[first3Norm==1]),
                    Movers = as.data.frame(datUse[pat_zip_1_0 != pat_zip_2_0]),
                    Age_18= as.data.frame(datUse[age06 <= 18]),
                    Same_Doc= as.data.frame(datUse[sameDocFirst2==1 & sameDocTwoThree ==1]),
                    Commercial = as.data.frame(datUse[commercial_1_0 == 1 & commercial_2_0 == 1 & commercial_3_0==1]),
                    MedicaidFFS = as.data.frame(datUse[(ffs_1_0) == 1 & (ffs_2_0) == 1 & ffs_3_0==1]))
    return(datList)
}

honoreDatList <- assignDatFunc(honoreDatUse)

##Get Number of women in each data set
nrowHonore <- lapply(honoreDatList, nrow)
for (nm in names(nrowHonore)){
    if(nrowHonore[[nm]] < minEstSize)honoreDatList[[nm]] <- NULL
}


honoreResList <- list()
honoreRatioList <- list()
##Define Honore Estimation
honoreEstFuncList <- list(distAndSwitch = function(dat) glm(one~ secondHospMinFirstHosp3 + time_current_diff- 1, family = 'binomial', data = dat))

##Honore Estimation
message('Honore Estimation')
for (typ in names(honoreEstFuncList)){
    honoreResList[[typ]] <- lapply(honoreDatList, honoreEstFuncList[[typ]])
    honoreRatioList[[typ]] <- lapply(honoreResList[[typ]], function(x)deltaMethod(x,"secondHospMinFirstHosp3/time_current_diff"))
}


##Print and save selected results
for (datTyp in names(nrowHonore)){
    outDir <- resultsDir
    listForOutput <-lapply(honoreResList[names(honoreResList) %like% 'Switch'], function(dat) dat[[datTyp]])
    stargazer(listForOutput, out = paste0(outDir,'Switch/',fileTyp,'_',datTyp,'_', 'feReg','.txt'), type = 'text',title = paste0('Data = ',gsub('_','-',datTyp)))
    stargazer(listForOutput, out = paste0(outDir,'Switch/',fileTyp,'_',datTyp,'_', 'feReg_','.tex'), type = 'latex',title = paste0('Data = ',gsub('_','-',datTyp)))
}


fileConn <- file(paste0(paramFilePref, "feHonoreSwitch.txt"))
writeLines(as.character(coefficients(honoreResList$distAndSwitch$Main)['secondHospMinFirstHosp3']), fileConn)
close(fileConn)

fileConn <- file(paste0(paramFilePref, "feHonoreDist.txt"))
writeLines(as.character(coefficients(honoreResList$distAndSwitch$Main)['time_current_diff']), fileConn)
close(fileConn)


message('Saving File')
save(honoreResList,honoreRatioList,file=paste0(datDir,'honoreFE_results_minPpl_',fileTyp,'_',minZipHospPpl,'cutAge',ageCut06,'.RData', sep = ''))


